Pairing and Superfluid Properties of Dilute Fermion Gases at Unitarity 
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We study the system of a dilute gas of fermions in 3-dimensions, with attractive interactions tuned 
to the unitarity point, using the non-perturbative Restricted Path Integral Monte Carlo (R-PIMC) 
method. The pairing and superfluid properties of this system are calculated at finite temperature. 
The total energy at very low temperature from our results agrees closely with that of previous 
ground-state Quantum Monte Carlo calculations. We identify the temperature T* ~ 0.70eF below 
which pairing correlations develop, and estimate the critical temperature for the superfluid transition 
T c ~ 0.25eF from a finite size scaling analysis of the superfluid density. 

PACS numbers: 



The recent experiments, starting with generating a 
degenerate gas of cold atoms 0, the use of Feshbach 
resonance to tune the effective interaction between the 
fermions 0, the measurements of the gap to identify 
a pairing scale and the measurements of vortices to 
identify the superfluid state, have all given an impetus 
to the study of superfluidity in the BCS-BEC crossover 
regime. 

The nature of superconductivity or superfluidity in a 
many-particle system with an increasing pairing attrac- 
tion was shown [l|,[2j to interpolate smoothly between the 
BCS regime for weak attraction between the fermions, 
and the Bose-Einstein condensation (BEC) regime for 
strong coupling. In the weak coupling regime, the Cooper 
pair size is much larger than the interparticle spac- 
ing, and the simultaneous pairing and condensation of 
fermions is well described by the BCS theory. In the 
strong coupling regime, fermions form strongly bound 
bosonic molecules at a pairing temperature scale T* that 
is considerably higher than the BEC temperature T c . 

In a two species system of fermions with an attractive 
interaction between them, the unitary point is defined 
by the divergence of the zero energy s-channel scattering 
length, a s , for two particles in free space. The unitary 
point is interesting since it describes a strongly correlated 
system with universal properties^. In general, there are 
two length scales that define the system: the interparti- 
cle spacing oc nT 1 / 3, , and the scattering length a s which 
contains information about the interaction potential. At 
the unitary point however, since a s diverges, all the prop- 
erties of the system are described by a single length scale 
kp 1 or energy scale 6f where kp = (Sir 2 ?t'n) 1 l z is the 
Fermi momentum of the corresponding non-interacting 
system. 



At the unitarity point there is no small parameter, 
therefore well controlled numerical methods are needed 
to calculate the properties of the system. Here we present 
the first calculation of the superfluid density and other 
pairing correlations as a function of temperature at the 
unitary point in a continuum model of a two compo- 
nent fermion gas with attractive pairwise interaction be- 
tween the species. Our main results are: (i) an accu- 
rate determination of the temperature dependent inter- 
nal energy which makes it a viable tool for thermometry 
for cold atoms; (ii) determination of the pairing scale 
T* ~ 0.7 'ep from growth of density correlations of oppo- 
site spin fermions; (hi) determination of the condensa- 
tion scale T c ~ 0.25ef from a finite size scaling analysis 
of the superfluid density. We use the restricted path in- 
tergal Monte Carlo (R-PIMC) technique, which is the 
fixed-node extension of the continuum PIMC method to 
fermionic systems [j| [ljj- This technique has recently 
been used to study helium-3, electron-hole liquids, and 
many body hydro gen [ll]. 

We consider an unpolarized system of two spins (or 
two hyperfine species) of particle density n. For the bare 
two-particle interaction at unitarity, we use the potential 
(also used by Q,) 
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where 2//i is the effective range of the potential. This 
potential has its only bound state at zero energy, with 
the eigenfunction given by tanh(/xr)/r, and a divergent 
scattering length a s . We measure distances in units of ro, 
the average interparticle spacing, which gives the density 
n = 3/(47rrg) ~ 0.238 and k F r = 1.919. The energy 
scale is eg = h 2 jmr\ and the Fermi energy ep — 1.841eo- 
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All the following results were obtained with pro = 12. 
In the dilute limit, with the interparticle spacing much 
larger than the range of the potential pr ^> 1, the exact 
form of the potential is unimportant; only the scattering 
length matters. 

Method: For a system of distinguishable quantum par- 
ticles, the density matrix in configuration space can be 
written as a path integral over coordinate space variables 
at discrete imaginary time intervals as, 



r r M 

p(Ro, R m ; /3)= dRi ■ ■ ■ dR M -i exp - S 



(2) 



Here R m denotes a configuration of N/2 up (labeled with 
t) and N/2 down (labeled with J,) fermions at the m th 
timeslice. The total extent in the imaginary time direc- 
tion, j3 = 1/ksT, is divided into M timeslices each of 
size r so that j3 — Mr. The link action S m is the sum of 
the kinetic-action terms and potential-action U m , given 
by, 



R*n 



2 4e r 
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(3) 



Eq. is the basis for a quantum-classical isomorphism 
between a system of quantum particles and a classical 
system of interacting polymers. Each particle path is 
interpreted as a polymer, with the beads of the polymer 
connected with springs described by the kinetic action. 
The potential-action describes the interaction between 
beads of different polymers. 

For indistinguishable fermions, the density matrix is 
anti-symmetrized over particle permutations, 



p F (R, R'; P) = -L^-iyppR, R>- (3) . 



(4) 



In the classical picture, a permutation V involving n par- 
ticles corresponds to the cutting and combining of the 
corresponding n polymers into one large polymer. As 
paths extend in space at low temperatures these quan- 
tum exchanges become more likely. 

For fermions however, a straightforward evaluation of 
the fermion density matrix given by Q leads to the 
fcrmion sign problem. It arises due to the cancellation, 
at low temperature, of approximately equal contributions 
from the positive and negative sign permutations, lead- 
ing to an exponentially vanishing signal-to-noise ratio in 
the Monte Carlo calculation 

The density matrix is a function of 2dN + 1 variables 
(the end configurations and time,) and can be solved for 
in a region of 'space-time' with restricted paths if the ini- 



tial condition pf(R, R'\ 0) 



J2 V S(VR- R'), and 



Dirichlet boundary conditions for r ^ 0, are fully speci- 
fied [Toj . It is natural to choose zero boundary conditions 
to define the path restriction region, so that knowledge 
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FIG. 1: The energy per particle in units of €fg = (3/5)eF as 
a function of T for N — 10 and 20 (dashed line; error bars not 
shown but comparable to N — 10) computed using the BCS- 
like trial density matrix ('BCS'). For comparison the energy 
obtained by using the non-interacting particle trial density 
matrix ('Free') is also shown. 



of the nodal surfaces of the density matrix allows com- 
putation of the thermodynamic properties of the system. 

To enable the computation, we utilize a trial den- 
sity matrix constructed to have reasonable physical and 
topological properties, that defines the nodal surfaces 
to restrict the paths. The R-PIMC method evades the 
sign problem for calculations which involve the diago- 
nal elements of the fermion density matrix. Due to 
the positivity of the density matrix on the diagonal, 
Pf(R, R;t) > 0, we can disregard any paths that start 
at a negative permutation of R since such paths have to 
cross a node at least once to end at R. We sum over 
only the positive (even) permutations, and keep only the 
positive node avoiding paths thus overcoming the sign 
problem. The permutation sum as well as the integrals 
in J3J are evaluated with the Metropolis Monte Carlo 
method. 

In the high temperature limit, it can be shown [l^ 
that the interacting density matrix is well approximated 
by the free particle density matrix, pp(R, R*;t) = 
(47re r)- dAr / 2 det[exp{-(r l -r*) 2 /4e r}]. In the low- 
temperature limit, the contribution of the excited states 
in the spectral expansion of the fermion density ma- 
trix, Pf {R,R';t) = exp(-r^„), are 

exponentially damped relative to the ground state. The 
density matrix can then be approximated for a non- 
degenerate ground state as, pf{R,R']t —> oo) = 
v I'o(-R) v I / o(^')- For the ground state trial wave function, 
we use a BCS-like antisymmetrized product of pairing 
functions |l3j . <fi(r^ — H-) between opposite spin fermions, 
and given by ^o(R) — det[r/>(r] —77)]. 



Results: The total energy vs. 
in Fig. Q]for N = 10 and N = 



temperature is shown 
20. The BCS-like trial 
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density matrix, in addition to its enabling the physics 
of pairing and condensation, is further justified because 
it gives the lowest total energy at low temperature of 
all the trial density matrices that we have tested. We 
also find good agreement of our energy extrapolated to 
T = with the value of E a /N = 0.44(l)e FG = 0.26e F 
obtained from the fixed node Green function Monte Carlo 
method[llEi, where tpo = 3/5ef is the average energy 
per particle of the non- interacting Fermi gas. 

Estimation of T* : The first indication of the crossover 
pairing scale T* « 0.70e F is obtained from a comparison 
of the energy obtained by using the free particle trial 
density matrix, and the BCS-like density matrix as shown 
in Fig. ^ For T > T* the free particle density matrix 
gives the lowest energy, whereas in the opposite regime 
T <T* the BCS density matrix gives the lowest energy. 

Further evidence of the pairing scale comes from g-\~[ 
pairing correlations as seen in Fig. [5] As the temperature 
is lowered, there is a strong enhancement of the on-site 
density correlations for opposite spins at T ~ 0.67eF (see 
inset) which identifies a pairing scale below which strong 
pairing correlations exist at temperatures well above the 
actual superfluid transition. 

In previous work[l7| a functional integral approach was 
used to estimate T* as a function of l/(fc F a s ) from a 
saddle point analysis of the gap and number equations. 
At unitarity they found T* ~ 0.57ei?- Upon including 
fluctuations to quadratic order around the saddle point, 
the transition was suppressed, especially around unitarity 
and beyond, to a lower T c . They estimated T c ~ 0.22e F 
at unitarity. The inclusion of the fluctuations around the 
saddle point to fourth order showed that the variations 
of T c is controlled by the coherence length kp^Q. In the 
BCS regime kp£o is exponentially large, while in the BEC 
regime, for a dilute Bose gas, kp^Q grows as a power law. 
In both these regimes, T c variations are small. However, 
around the unitary point, fc F £o ~ and the fluctua- 

tions in T c /ep are of order unity. This is precisely where 
simulations are of greatest value since the system is in 
a strongly correlated regime inaccessible to perturbation 
theory. 

We compute the superfluid density to directly esti- 
mate the critical temperature T c for the superfluid transi- 
tion. The superfluid density p s for an A-particle system, 
shown Fig. |3J is calculated from the winding number 
estimator 19], 



p 



(W 2 ) 

2e (3N 



(5) 



Here, W is the winding number defined as the number of 
times periodic boundary conditions are invoked as paths 
start from Ri and end at a periodic image R-pi. 

We start by reviewing the Josephson relation which 
is a hyperscaling relation between the correlation length 
exponent v defined by £ ~ t~ u and the superfluid density 
exponent £ defined by p s ~ ft, where t — (T — T c )/T c 




FIG. 2: The pair correlation function of the two species of 
fermions at different values T. Note the sharp increase in pair 
correlation from 0.67T F to 0.54T F . 
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FIG. 3: The superfluid fraction p s /p as a function of T for 
N = 20. 



is the reduced temperature. The singular part of the 
free energy density near the transition is f s (t) ~ £,~ d , 
and since p s (V(/>) 2 ~ f s , we obtain p s ~ £-( d ~ 2 ) ^ 
iv(d-i) which directly gives the Josephson relation 2CI 
Q/v = d—2. In a finite system, the scaling hypothesis 2fJ 
f s (t,L) = L~ d T[L/i{t)] states that the free energy de- 
pends on t only through the ratio of the system size L 
and the bulk correlation length £. Using the Josephson 
relation, this implies that p s (t,L) = L 2 ~ d Q[L/£(t)]. In 
d = 3, linearizing the function Q near t = and using 
L - iV 1 / 3 leads to 



Q 



L 



Ps(t) 



Q(0) + qN 1/3v 



T-T c 
T c 



(6) 



where q is a constant. In Fig. 0]we plot N 1 / 3 p s (t)/ p vs T 
for several system sizes N. At the transition temperature 
T c , the size dependence vanishes and all the curves meet 



4 







N=10 
N=20 
N=30 




- 












X 






X 
















- 






- 








+ ~~ — - 


1 1 


1 1 1 


i i 


i'* 



0.15 0.175 0.2 0.225 0.25 0.275 0.3 0.325 
T[E F ] 

FIG. 4: The scaled, linearized, superfluid density, N 1 ^ 3 p a /p, 
as a function of T for N = 10, 20, and 30. The intersection of 
the three lines gives the superfluid transition temperature as, 
T c « 0.25T F . 

at a point which determines the critical temperature T c ss 
0.25e F . Our result agrees well with the estimate of T c w 
0.22 ep 01 obtained by including fluctuations around the 
saddle point. 

The transition temperature has also been calculated 
using lattice Monte Carlo techniques 0, H3|. Our es- 
timate is however higher than the lattice Monte Carlo 
estimate of Burovski et al. of T c w 0.15 ef 22J . The 
main distinction of our path integral Monte Carlo is that 
we work directly in the continuum so the unitary limit is 
perfectly well-defined. The lattice simulations have to ex- 
trapolate T c to the zero filling factor limit in order to get 
the correct behavior in the unitary limit. On the other 
hand, for the attractive-U Hubbard model there is no 
sign problem|l6j which is a definite advantage. However, 
given the good agreement of our finite temperature R- 
PIMC method with the zero temperature GFMC results 
for the energy, we believe we have an accurate description 
of the nodal surface. The R-PIMC method can further 
be used to study the competition between superfluidity 
and magnetization in unequal fermion populations. 

One of the difficulties in the experiments is to de- 
termine the temperature precisely. Usually a Maxwell 
Boltzmann fit to the excited atoms yields an estimate. 
However, at low temperatures when the number of atoms 
in the excited states is greatly reduced such an estimate 
becomes unreliable ^^|. The strong dependence of the en- 
ergy on temperature above T c seen from our results, indi- 
cates that measurements of the mean field energy can be 
converted to a temperature scale, with appropriate cor- 
rections for a trap using a local density approximation. 
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